function [hyp1cov] = gptrain(hyp,covfunc,x,y,params)
%GPTRAIN	optimise hyperparameters based on the negative log marginal likelihood
%	[hyp1] = gptrain(hyp,covfunc,x,y1,params)
%	
%
%	Inputs:     hyp         structure with lik and covf hyperparameters
%               covfunc     covariance function to be used
%               x, y1       input, output variables
%               params      parameters to normalise (e.g.:
%                           params = [1 0 5 .5; ...  % which parameter to optimise; lower-end of interval; upper-end of interval; step within interval;
%
%	Outputs:   hyp         covf hyperparameters selected
%		
%		
%	
%	See also demogpsimilarity

%	References: Gaussian Process for Machine Learning Book.
%	
%	

%	Copyright 2013 MAF Pimentel
%	This program is free software: you can redistribute it and/or modify
%	it under the terms of the GNU General Public License as published by
%	the Free Software Foundation, either version 3 of the License, or
%	(at your option) any later version.
%	
%	This program is distributed in the hope that it will be useful,
%	but WITHOUT ANY WARRANTY; without even the implied warranty of
%	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%	GNU General Public License for more details.
%	
%	You should have received a copy of the GNU General Public License
%	along with this program.  If not, see <http://www.gnu.org/licenses/>.


%	$Author: MAF Pimentel$
%	$Revision: 1.0.0.0$
%	$Date: 11-Jun-2013 11:35:49$
%	Contact: marco.and.pimentel@gmail.com
%	Originally written on: PCWIN64

N = size(params,1);
%fprintf('\nOptimising %d parameters\n',N);

for n = 1 : N
    p{n} = params(n,2):params(n,4):params(n,3);
    p{n} = [params(n,1) p{n}];
end

switch N
    case 1
        hyp1cov = hyp.cov;
        nlmlK = NaN*ones(length(p{1})-1,1);
        for i = 1 : length(p{1})-1
            hyp.cov(p{1}(1)) = p{1}(i+1); 
            [hyp, ~, ~] = ...
                minimize(hyp, @gp, -200, @infExact, [], covfunc, @likGauss, x, y);
            nlmlK(i) = gp(hyp, @infExact, [], covfunc, @likGauss, x, y);
        end
        
        [~,c] = min(nlmlK);
        hyp1cov(p{1}(1)) = mean(p{1}(c+1));
    
    case 2
        hyp1cov = hyp.cov;
        nlmlK = NaN*ones(length(p{1})-1,length(p{2})-1);
        for i = 1 : length(p{1})-1
            hyp.cov(p{1}(1)) = p{1}(i+1);
            for j = 1 : length(p{2})-1
                hyp.cov(p{2}(1)) = p{2}(j+1); 
                [hyp, ~, ~] = ...
                    minimize(hyp, @gp, -200, @infExact, [], covfunc, @likGauss, x, y);
                nlmlK(i,j) = gp(hyp, @infExact, [], covfunc, @likGauss, x, y);
            end
        end
        
        [r,c] = find(nlmlK == min(nlmlK(:)));
        hyp1cov(p{1}(1)) = mean(p{1}(r+1));
        hyp1cov(p{2}(1)) = mean(p{2}(c+1));
        disp('');
        
    otherwise
        fprintf('\nCode is not ready to optimize more than two variables\n');
        return;
end

end